/*
Output: Figure 2. Average Annual Percent Area Harvested, Enrolled and All Parcels
*/


run "/Users/holyfantastic/Dropbox/Forest/Paper/Nature/FinalDocuments/Code/000declare_path" // change to your local path

set scheme white_jet

use "ParcelLevelDataSept282023_VerNov82024_R1.dta" , clear

replace Landtr_Cut_mean_forest = Landtr_Cut_mean_forest *100

replace YearsOwned =. if YearsOwned < 0
replace dis_road = dis_road * 1000
label var YearsOwned ""

preserve
	keep if group == 3
	gen group_relabel = 1
	save "$temp/temp" , replace
restore

gen group_relabel = 2
label define lab_group_relabel 1 "Enrolled Parcels" 2 "All Parcels"
label values group_relabel  lab_group_relabel

append using "$temp/temp"

******************************************************************************

local i = 1

foreach y in  Landtr_Cut_mean_forest8599 Landtr_Cut_mean_forest0014 Landtr_Cut_mean_forest1519{
	
	preserve
	
		collapse (mean) mean_`i'= `y' (sem) sem_`i'= `y', by(group_relabel)
		gen lower_bound_`i' = mean_`i' - 1.96 * sem_`i'
		gen upper_bound_`i' = mean_`i' + 1.96 * sem_`i'

// 		collapse (mean) mean_`i'= `y' (sd) sd_`i'= `y', by(group_relabel)
// 		gen lower_bound_`i' = mean_`i' - 1 * sd_`i'
// 		gen upper_bound_`i' = mean_`i' + 1 * sd_`i'
		
		gen x = `i'
		save "$temp/temp_`i'" , replace
		
	restore
	
	local i = `i' + 1

}

clear

append using "$temp/temp_1"
append using "$temp/temp_2"
append using "$temp/temp_3"

gen x2 = _n
replace x2 = x2 - 0.5 if group_relabel == 2

tw || ///
	(bar mean_1 x2 if group_relabel == 1, barwidth(0.5) fcolor(gs8) lcolor(gs8) lwidth(medium) legend(off)) ///
	(bar mean_2 x2 if group_relabel == 1, barwidth(0.5) fcolor(gs8) lcolor(gs8) lwidth(medium) legend(off)) ///
	(bar mean_3 x2 if group_relabel == 1, barwidth(0.5) fcolor(gs8) lcolor(gs8) lwidth(medium)) ///
	(rcap lower_bound_1 upper_bound_1 x2 if group_relabel == 1, lw(medium) lcolor(black) lp(solid)) ///
	(rcap lower_bound_2 upper_bound_2 x2 if group_relabel == 1, lw(medium) lcolor(black) lp(solid)) ///
	(rcap lower_bound_3 upper_bound_3 x2 if group_relabel == 1, lw(medium) lcolor(black) lp(solid)) ///
	(bar mean_1 x2 if group_relabel == 2, barwidth(0.5) fcolor(gs12%30) lcolor(gs8) lwidth(medium) legend(off)) ///
	(bar mean_2 x2 if group_relabel == 2, barwidth(0.5) fcolor(gs12%30) lcolor(gs8) lwidth(medium) legend(off)) ///
	(bar mean_3 x2 if group_relabel == 2, barwidth(0.5) fcolor(gs12%30) lcolor(gs8) lwidth(medium)) ///
	(rcap lower_bound_1 upper_bound_1 x2 if group_relabel == 2, lw(medium) lcolor(black) lp(solid)) ///
	(rcap lower_bound_2 upper_bound_2 x2 if group_relabel == 2, lw(medium) lcolor(black) lp(solid)) ///
	(rcap lower_bound_3 upper_bound_3 x2 if group_relabel == 2, lw(medium) lcolor(black) lp(solid)) ///	
	,  note("") legend(on order( 1  "Enrolled Parcels"  7 "All Parcels" ) pos(1) rows(1))  xtitle("")  ///
	ylabel(1(0.5)3 , nogrid) ytitle("Percent Area Harvested") xlabel(  1.25 "1985-1999" 3.25 "2000-2014" 5.25 "2015-2019"  , labsize(small) angle(45) nogrid) 

	foreach fig_path in figure {
	
	 dis  "$`fig_path'/fig2_bar_four_average_harvest_R1"
	 
	graph export "$`fig_path'/fig2_bar_four_average_harvest_R1.png", as(png) name("Graph") replace
	graph export "$`fig_path'/fig2_bar_four_average_harvest_R1.pdf", as(pdf) name("Graph") replace

}
